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■ Three numerical algorithms are proposed to solve the time-dependent elastodynamic equations in 

elastic solids. All algorithms are based on approximating the solution of the equations, which can be 
written as a matrix exponential. By approximating the matrix exponential with a product formula, 
an unconditionally stable algorithm is derived that conserves the total elastic energy density. By 
expanding the matrix exponential in Chebyshev polynomials for a specific time instance, a so-called 
"one-step" algorithm is constructed that is very accurate with respect to the time integration. By 
formulating the conventional velocity-stress finite-difference time-domain algorithm (VS-FDTD) in 
matrix exponential form, the staggered-in-time nature can be removed by a small modification, and 
higher order in time algorithms can be easily derived. For two different seismic events the accuracy 
of the algorithms is studied and compared with the result obtained by using the conventional VS- 
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An important aid in the understanding of wave propagation in inhomogeneous media is seismic forward modeling. 
CIh' In all but the simplest cases, an analytical solution of the elastodynamic equations is not available, and one must 
resort to numerical solutions. For this, two main strategies can be followed: one solves the elastodynamic equations 
either in the strong formulation, where the equations of motion and boundary conditions are written in differential 
form, or in the weak formulation, where the equations of motion are given in integral form. The latter formulation, 
implemented by finite element [1,2], spectral element [3,4] or finite integral methods [5,6], may be preferred to deal 
, with complex geometries, or non-trivial free surface boundary conditions. 

The first strategy (the strong formulation) is followed in the finite-difference approach that is based on solving either 
the first-order velocity-stress differential equations [7-9], or the second order wave equation [10-12]. In the original 
formulation of the velocity-stress finite-difference time-domain (VS-FDTD) approach [7,8], space and time are both 
discretized using second-order finite differences. Many enhancements have been introduced to increase the accuracy 
. and treatment of special boundary conditions. For example, spectral methods have been employed to increase the 
q ' accuracy of the approximation of the spatial derivative operators [13,14]; a rotated staggered spatial grid has been 
developed to surmount some instability and lack of spatial accuracy issues [15]; polynomial expansions of the time 
evolution operator have been used to increase the accuracy and efficiency of the time integration process [14,16]. 
Although each specific method cited above offers its own advantages and drawbacks, none of these methods feature 
Qh unconditional stability and/or exact energy conservation, which can be useful for long integration times, and high 
material contrast situations, where instabilities have been noticed [17]. The algorithms that will be introduced in this 
paper, address this issue. 

Starting from the velocity-stress first-order differential equations, it can be shown (see section II) that the solution of 
these equations can be written as the matrix exponential of a skew-symmetric matrix. This constitutes an orthogonal 
transformation, conserving the total energy density. Guided by recent results regarding the numerical solution of 
the time-dependent Maxwell Equations [18-23], where the underlying skew-symmetry of the equations of motion 
is exploited in a matrix exponential approach, we apply this framework to the current problem. In section III 
this framework is briefly repeated for convenience and three algorithms are derived to solve the time-dependent 
elastodynamic equations. The incorporation of the presence of a source is described in section IV. For some typical 
examples, the performance and efficiency of the algorithms is studied in section V, and the conclusions are summarized 
in section VI. 



1 



II. THEORY 



In the absence of body forces, the linearized equation of momentum conservation reads [24] 

j J 

where p is the density, Ui is the displacement field and <jij the stress field (i = x, y, z). It can be recast into a coupled 
first order velocity-stress equation [7] , yielding in matrix form 



8_ ( a \ _ ( -CD T \( a 
dt { v ) ~ UD v 



(2) 



Here, a = {a xx ,u X y,u xz ,a vx ,Uyy,Uy Z ,a zx ,a zy ,a zz ) T ), v = (v x ,v y ,v z ) T is the velocity field and D is the matrix 
containing the spatial derivatives operators, 

dx 2 dy 2 dz 2 dy U 2 dz U U \ 

D= \l \l £ |£ ±£ . (3) 

\ u u 2 dx u ^ 2 dy 2 dx 2 dy dz J 

The stiffness tensor C = {Cijki} relates the stress and strain, 

°rj = CijklCkh (4) 

and is symmetric and positive definite for elastic solids. 

Some important symmetries in matrix equation (2) can be made explicit by introducing the fields 

w = Vpv, (5) 

and 

s = J=.. (6) 

The expression \fC is valid since C is symmetric and positive definite. 
By definition, the length of tjj = (s,w) T , given by 

II^H 2 = (V#) = / iJj T ^dr= f (w 2 + s 2 )rfr= f (pv 2 +<7 T C- 1 <j)dv (7) 
Jv Jv Jv 

is related to the elastic energy density 

w = \ (pv 2 + o T C~ x a) = i j pv 2 + ^2 Cijkieijeki j , (8) 

of the fields. 

In terms of ip, matrix equation (2) becomes 



d ( -VCD T ±\ 
Using the symmetric properties of p and \[C \ one can prove that the matrix H is skew-symmetric 

/ o (-tsWct) \ / o Vcd 



H 1 



V 



with respect to the inner product as defined in equation (7). The formal solution of equation (9) is given by 

^(t) = e t? V(0)=W(tM0), (11) 



2 





















— ^ 


V J< [ 

* 


CT ZX 




FIG. 1. Unit cell of the three-dimensional staggered grid onto which the continuous velocity and stress fields of the elas- 
todynamic equations are mapped in order to conserve the skew-symmetry. Left: grid for elastic isotropic solids. Note: the 
Lame constants A and /i coincide with the stress field components, and the mass density is only defined on velocity field points. 
Right: grid for the general anisotropic case (i,j = x,y,z). 



where ^(0) represents the initial state of the fields and the operator U determines their time evolution. And, since 7i 
is skew-symmetric the time evolution operator, U is an orthogonal transformation: 

U(t) T = U{-t) = U-\t) = e~ m , (12) 

and it follows that 

(u(t)m\mm) = (mm)) = mmm- as) 

Hence, the time evolution operator U(t) rotates the vector ip(t) without changing its length \\ip\\- In physical terms, 
this means that the total energy density of the fields does not change with time, as can be expected on physical 
grounds [24]. 

In practice, the construction of a numerical algorithm requires to discretize space and time. During both these 
procedures, the skew-symmetry of TL (during the spatial discretization) and the orthogonality of U (during the time 
integration) should be conserved. For the discretization of space, this requirement can be met by choosing a staggered 
spatial grid [8] and a central difference approximation for the spatial derivative. This yields a skew-symmetric matrix 
H for the discrete analogue of TL. A unit cell of the grid is shown in figure 1. The explicit form of H is derived in the 
appendix, in the case of a two-dimensional isotropic elastic solid. Accordingly, the discrete analogue of ij)(t) is given 
by vector *£(t). 

The continuous problem, defined by TL, is now translated to a lattice problem defined by H: 

= exp(i/f)^(0) = J7(t)*(0). (14) 

or, in the time-stepping approach, we have for a small timestep r 

*(* + t) = exp(rfl")*(t) = U(r)^(t). (15) 

At this point, we invoke three different strategies to perform the time integration, i.e. to approximate the matrix 
exponential exp(iiJ). Here, we closely follow the derivation of algorithms to solve the time-dependent Maxwell 
equations [18-23], where the problem to be solved is stated in a very similar form, although the underlying physics is 
different. The first algorithm is based on conserving the existing symmetries during the discretization of time, and is 
unconditionally stable. Here, the time integration is carried out by a time-stepping procedure. The second algorithm 
is based on approximating the solution itself for a particular time instance, by means of a Chebyshev expansion, 
and constitutes therefore a "one-step" algorithm. The last algorithm is based on recasting the original velocity-stress 
finite-difference algorithm into matrix exponential form. This allows to remove the staggered-in-time nature, and 
offers an elegant way to derive higher-order in time algorithms. The construction of the algorithms is briefly repeated 
in the next section. 

III. THE MATRIX EXPONENTIAL APPROACH 
A. Unconditionally stable algorithms 

A sufficient condition for an algorithm to be unconditionally stable is that [25] 

\\u(T)*(t)\\ < (i6) 
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Since U{t) is an orthogonal transformation (see previous section), we have ||[/(t)\E'(£)|| = ||\&(t)||, and it is sufficient to 
conserve the orthogonality of U{t) for an approximation U(t) to U(t), in order to construct an unconditionally stable 
algorithm. One way to accomplish this is to make use of the Lie- Trotter-Suzuki formula [26,27] and generalizations 
thereof [28,29]. If the matrix H is decomposed, so that H — YTi=i then 

f/i(r) = e rHl ...e rH ", (17) 

is a first order approximation to U(t). More importantly, if each matrix Hi is skew-symmetric, then U\(t) is orthogonal 
by construction, and hence, algorithms based on U\{t), are unconditionally stable. Using the fact that both U(t) 
and U\{t) are orthogonal matrices, the error on U\{t) is subject to the upper bound [30] 

||^T)-lMr)||<y£||[i^]||, (18) 

i<j 

where [H h Hj] = HiHj - H 3 H t . 

In the appendix, the decomposition of H is carried out for two-dimensional isotropic elastic solids, for which p = 12, 
and it is shown that each matrix H t is block diagonal. The computation of the matrix exponential of a block diagonal 
matrix Hi can be performed efficiently, as it is equal to the block-diagonal matrix of the matrix exponentials of the 
individual blocks. Therefore, the numerical calculation of e rHi reduces to the calculation of matrix exponentials of 
2x2 matrices, which are rotations. 

In practice, implementation of the first order algorithm is all that is required to construct higher order algorithms. 
This is due to the fact that in the product-formula approach, the accuracy of an approximation can be improved in a 
systematic way by reusing lower order approximations, without changing the fundamental symmetries. For example, 
the orthogonal matrix 

U 2 (t) = U 1 (-t/2) T U 1 (t/2) =e TH "/ 2 ...e rHl / 2 e rHl / 2 ...e TH "/ 2 , (19) 

is a second-order approximation to U{t) [28,29]. A particularly useful fourth-order approximation (applied in for 
example [18,19,27-39]) is given by [28] 

U 4 (t) = U 2 (aT)U 2 (aT)U 2 ((l - 4a)r)C/ 2 (ar)C/ 2 (ar), (20) 

where a = 1/(4 -4 1 / 3 ). 



B. One-step algorithm 



A well-known alternative for timestepping is to use Chebyshev polynomials to construct approximations to time- 
evolution operators [40-44] . This approach has also been successfully applied to the problem of seismic wave propaga- 
tion [14,16]. However, the main differences between these implementations and the algorithm explained below, besides 
the spatial discretization (which is based here on a central difference approximation, instead of a spectral method), 
is that in the present case, matrix H is explicitly anti-symmetrized, which gives rise to purely imaginary eigenvalues 
[45]. This is an important property to justify the validity of the expansion. In the derivation of the current algorithm, 
we follow [20-22], a recent implementation of the Chebyshev algorithm to solve electromagnetic wave propagation. 

The basic idea is to expand the time evolution matrix U (t) = exp(iTJ) for a specific time instance t in matrix valued 
Chebyshev polynomials on the domain of eigenvalues of H, which lies entirely on the imaginary axis since H is skew- 
symmetric. For proper application of the expansion, the domain of eigenvalues is rescaled to [—1,1], by considering 
the matrix B = —iH/\\H\\i, where ||i?||i denotes the 1-norm of the matrix. It is given by \\H\\i = max., J2i \^ij\, see 
[45], and is easy to compute since the matrix H is sparse. Operating on state \& (0) , the expansion becomes 



tf (t) = cxp(ti7)*(0) = exp(izB)*(0) 



Mz)I + 2^r t J n (z)T n (B) 



1=1 



*(0), 



(21) 



where / is the identity matrix, z = t\\H\\, J n are nth order Bessel functions and T n {B) 
Chebyshev polynomials, defined by the recursion relation 



i n T n (B) are the modified 



T o (B)*(0) = tf(O), 
Ti(B)*(0) =iB*(0), 
f„ +1 (B)*(0) = 2iST n (B)*(0) +f„_i(B)*(0) for n > 1. 



(22a) 
(22b) 
(22c) 
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Due to the fact that the matrix B is purely imaginary, it follows from the above recursion relation (22) that T n (B)^> (0) 
and thus ^(t) will be real valued and no complex arithmetic is involved, as should be the case. 

In practice, the summation in Eq. (21) will be truncated at some expansion index m. This number depends on 
the value of z, since the amplitude of the coefficients J n {z) decrease exponentially for n > z; this is explained in 
more detail in Refs. [20-22]. Consequently, the computation of one timestep amounts to carrying out m repetitions 
of recursion relation Eq. (22) to obtain the final state. This is a simple procedure: only the multiplication of a vector 
with a sparse matrix and the summation of vectors are involved. 



C. A modified VS-FDTD algorithm 

In section III A, it was shown that in the product formula formalism, higher-order in time algorithms can be 
constructed by reusing the lower order algorithms. This elegant technique to increase the accuracy of time integration 
can also be applied to FDTD algorithms (in case of the Maxwell Equations, see [23]), and hence the conventional 
VS-FDTD algorithm, if it is recast into an exponent operator form. 

The update equations of the VS-FDTD [7,8] algorithm can be written as 

: { ?: T % ) - e + + rm ( 'w ) , u M ( ) , (23) 



where 



and 



-CD T 




(24) 



B =(i!)' (25) 

Since A 2 = B 2 = 0, the time evolution operator U\{t) is equal to 

Ui{t) = cxp(rA) exp(rB), (26) 

and second-order accurate in time operating on fields that are defined staggered- in-time. However, if U\(t) is inter- 
preted as an approximation to the operator exp(-7vl + tB), working on fields non staggered-in-time, it is a first-order 
approximation. Furthermore, since the time evolution operator is now expressed in matrix exponential form, the 
accuracy can be increased by the same procedure as was used for the unconditionally stable algorithms (cf. eqs. (19) 
and (20)). Therefore, the operator 

U 2 (t) = exp(rA/2) exp(rB) exp(r,4/2), (27) 

constitutes a second-order approximation to the exact time evolution operator. And similarly, a fourth order algorithm 
can be derived, see Eq. (20). 

So, by introducing a small modification to the original VS-FDTD algorithm, that would require a minimal change 
in existing numerical codes, the staggered-in-time nature is removed, and higher-order in time algorithms are derived. 

IV. SOURCES 

In the presence of an explosive initial condition, or other time-dependent body force, the equation of motion reads 

^(t) = nm + <P(t), (28) 

where the time-dependent source is denoted by the term <j>(t). The formal solution is given by 

tp(t) = exp(tH)ip{0) + [ exp((i - u)H)<j>(u)du. (29) 
Jo 
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In the time-stepping approach (e.g. the unconditionally stable algorithms), the source term -if its time dependence is 
known explicitly- can be integrated for each timestep. For example, a standard quadrature formula can be employed 
to compute the integral over u, like the fourth-order accurate Simpson rule [47] 

J' +T e {t+T -^ n ^{u) du^^ (e™ cj>(t) + Ae TH ' 2 4>(t + r/2) + <j>{t + r)) . (30) 

In case of the one-step algorithm, this approach of incorporating the source term is not efficient, as for each value 
of t — u the recursion (22) would have to be performed, and a different route is taken. Instead, each of the two terms 
on the right hand side of equation (29) is expanded in modified Chebyshev polynomials separately. The expansion of 
the first term is discussed in section IIIB. For the second term, the source-term integral, the procedure is carried out 
here for a source with Gaussian time dependence, 

0(t) = /(t)S(r) = exp(-a(t - t ) 2 )S(r), (31) 

where S(r) denotes the spatial dependence of the source. 
The source term 

h(t,t ,a,H)= f due {t ~ u)H '/(it), (32) 
Jo 

is expanded in modified Chebyshev polynomials, 

h(t, t , a, H)S(r) = (^b I + ^ b k f^j S(r), (33) 

where the expansion coefficients are given by 
2 f* 

b k = i~ k - h(t, t , a, cos 9) cos k9d9. (34) 
Jo 

The replacement of H by cos 9 emphasizes that H should be normalized such that all eigenvalues lie in the range 
[—1,1]. We proceed by evaluating h(t,to,a,H). After substitution of z = ||-ff||i, zo — to||.ff||i, ft — a/\\H\\f and 
x = —iH/\\H\\i normalize the matrix H, we obtain 



h(z,z ,@,x) = ^y^ ex P (( z - z o) ix - x 2 /4f3) erf (^zo\/~P - + erf (j z ~ z o)v^ + 



IX 



(35) 



2^y. 

Now we put x — cos 9 and the remaining integral over 9 in equation (34) is computed by a Fast Fourier transformation: 
b k = 2 l - k e 27rmfc/Ar M^^o,/3,cos^). (36) 

n=0 

The derivation of the Chebyshev expansion coefficients for a source defined by equation 

g(t) = ^/(t) = -2a(t - t ) exp(-a(t - to) 2 ). (37) 
is very similar, and will not be treated explicitly here. 



V. RESULTS 



The performance and accuracy of the algorithms introduced in the previous sections is studied by comparing the 
results with a reference solution generated by the one-step algorithm, denoted by $(t). This choice is motivated by 
the fact that the latter, considering the time integration, produces numerically exact results [14,16,41]. Furthermore, 
there are rigorous bounds on the error of the unconditionally stable algorithm (cf. Eq. (18)): in the presence of a 
source <j>(t), the difference between the exact solution &(t) and the approximate solution &(t), obtained by using the 
4th order unconditionally stable algorithm, is bounded by [23] 

||*(t) - *(*)» < c 4 tr 4 (||*(0)|| + jT duU(u)\\^ , (38) 
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FIG. 2. The corner-edge system, consisting of two different materials. The system size and location of the source (x) are 
indicated in the picture. At the top, a free-surface boundary condition is imposed, the other boundaries are rigid. The overall 
density is p — 2500 kgm -3 , and the mesh size is S = 100 m. In the bulk material (I), the wave velocities are v p — 6 kms -1 and 
v s = 2 kms -1 , whereas in the inner material (II), they are v p — 9 kms -1 and v s = 3 kms" 1 . The source excites the s xx and 
s yy stress fields with time dependence g(t) from equation (37) and parameters a = 40 and to = 1.5 s. 




FIG. 3. Two snapshots of the kinetic-energy density distribution of the corner-edge system of figure 2. Left: state at t = 1.8 
s, right: state at t = 6.0 s. 



where C4 is a constant. For the difference between the exact solution and the solution obtained by using the one-step 
algorithm, we can write using the triangle inequality 

||*(f) - *(*)|| < \Mt) *(*)« + ll*(*) - (39) 

Using equations (38) and (39) and the fact that the difference 11^(4) — ^(4)|| vanishes with r 4 , as we will show below, 
we can be confident that the one-step algorithm indeed produces numerically exact results. This justifies to define 
the error as ||*(4) - #(4)||/||*(4)||). 

The performance of the following algorithms is considered: the original VS-FDTD algorithm, denoted by Vir; the 
unconditionally stable algorithms, denoted by LTS-2 and LTS-4, for resp. second and fourth order accuracy in time; 
the non staggered-in-time, modified VS-FDTD algorithm, denoted by VNS-2 and VNS-4, depending on the accuracy 
in time. 

Consider a rectangular system consisting of two different materials, displayed in figure 2. This system is also 
studied in references [8,10], and proved to be a good testing situation for the performance of an algorithm solving the 
elastodynamic equations. The time evolution of the velocity and stress fields is computed using an explosion as initial 
condition, modeled by equation (37), up to t = 6 s, with the six different algorithms. In figure 3, the kinetic-energy 
density distribution is shown for two different time instances. 

For the corner-edge system, the errors are listed in table I, and are shown in figure 4. For the largest timestep, 
the VS-FDTD algorithms (Vir,VNS-2,VNS-4) are unstable. This can be expected, since the maximum timestep is 
limited by the largest velocity, the v p velocity, and the mesh-size, through the Courant limit [8] 

T <^- ( 4 °) 

Furthermore, from table I it is clear that for all algorithms the error scales according to the order of accuracy in 
time. We also see that for the corner-edge system, the VS-FDTD algorithms perform much better than the energy- 
conserving LTS algorithms, as long as the timestep is smaller than the Courant limit. For timesteps larger than the 
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Time-step (s) 


Vir 


VNS-2 


VNS-4 


LTS-2 


LTS-4 


5 ■ 1(T 2 


oo 


oo 


oo 


3.0 ■ 10 u 


3.1 ■ 10 u 


5 ■ 10" 3 


1.7 ■ 10 -2 


1.7 ■ 1(T 2 


4.1 ■ 1(T 6 


5.9 ■ 10" 1 


1.6 ■ 10 -4 


5 ■ 10" 4 


1.7- 1CT 4 


1.7- 10~ 4 


4.1 • 10 -10 


6.1 ■ 10~ 3 


1.6 ■ 10 -8 


5 • 10" 5 


1.7 ■ 10" 6 


1.7 ■ 1(T 6 


2.0 ■ 10" 13 


6.1 ■ 10" 5 


2.0 ■ 10~ 10 


5- 10" 6 


1.7 ■ 1CT 8 


1.7 ■ 10~ 8 


2.8 ■ 10" 13 


6.1 ■ 10~ 7 





TABLE I. Error as function of timestep for all timestepping algorithms, for the system defined in figure 2. An infinite symbol 
(oo) denotes that the algorithm was not stable, whereas in one case (-) the computation was not performed. The error in the 
staggered-in-time Vir algorithm is determined by averaging the error in the kinetic and potential energy density at respectively 
the final time instance and the final time instance shifted by half a timestep. The time shifting procedure is carried out by the 
Chebyshev algorithm and also applied to prepare the initial condition. 
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i_ 1e+00 
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1e-04 
1e-06 



1e-08 



1e-10 



1e-12 



1e-14 

1e-06 1e-05 1e-04 1e-03 1e-02 1e-01 

timestep (s) 

FIG. 4. Error as function of timestep for all timestepping algorithms, for the system defined in figure 2, using the data from 
table I. 



Courant limit, the LTS algorithms (LTS-2 and LTS-4) are stable, although the error does not (yet) scale according 
to the order of accuracy in time [48]. For very accurate results (errors below 10 -10 ), the number of operations the 
achieve this accuracy becomes so large that the error does not scale systematic anymore. 

With respect to the efficiency of the algorithms, we note that the number of matrix- vector operations W, necessary 
to perform one timestep, is 1 for the Vir algorithm, 1.5 for the second-order VNS-2 and LTS-2 algorithms, and 10 
for the fourth-order VNS-4 and LTS-4 algorithms. For the specific example here, the corner-edge system, the one- 
step algorithm employs m = 1514 expansion terms. At t = 6 and r = 0.005, the VNS-2 algorithm already uses 
more (namely 1800) matrix-vector operations. Therefore, we draw the conclusion that one-step algorithm should be 
preferred to be used to solve the time evolution for this problem. Note that in general, the choice of which algorithm 
to use depends heavily on which degree of error is acceptable. In this specific case, there are values for r for which 
the VNS-2 algorithm uses less matrix-vector operations than the one-step algorithm, but then the error will be larger 
than the error for r = 0.005, and maybe unacceptably high. On the other hand, one VNS-2 or LTS-2 matrix-vector 
operation is carried out (in practice) faster than one Chebyshev recursion iteration, although this depends on the 
actual implementation. 

It is important to note that the initial condition plays an important role in the error of the solution produced by 
a specific algorithm. From the results of the corner-edge system, one might draw the conclusion that the VS-FDTD 
algorithms achieve better results than the LTS algorithms for all systems. This is not true. In table II, the error is 
listed as function of timestep for all algorithms, as compared with the one-step algorithm, for a system consisting of a 
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Time-step (s) 


Vir 


VNS-2 


VNS-4 


LTS-2 


LTS-4 


5 ■ 1(T 2 


oo 


oo 


oo 


1.2 ■ 10 u 


5.8 ■ 10-^ 


5 ■ 1(T 3 


1.3 ■ 10 _1 


i.3 ■ kt 1 


1.6 ■ 1(T 5 


5.4 ■ 10~ 2 


1.3 ■ 10 -5 


5 ■ 10" 4 


1.3 ■ 1CT 3 


1.3 ■ 1(T 3 


1.6 ■ 1(T 9 


5.4- 1(T 4 


1.4 ■ 1(T 9 


5 • 10~ 5 


1.3 • 1CT 5 


1.3 ■ 10~ 5 


1.6 ■ 10" 13 


5.4 ■ 10~ 6 


1.7- lO" 11 


5- 10" 6 


1.3 • 10" 7 


1.3 ■ 10~ 7 


1.1 ■ 1CT 13 


5.4 ■ 10™ 8 


1.7- 1CT 10 



TABLE II. Error as function of timestep for all timestepping algorithms. In this case, the one-step algorithm employs 
m — 321 expansion terms. The system measures L x = L y — 10, with a mesh of S — 0.1, and the material parameters p, A,/i 
vary randomly in space with values distributed randomly in the interval [1, 3]. The error is determined at t = 3. (all quantities 
are expressed here in dimensionless units). 
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1e-05 



1e-04 



1e-03 1e-02 

timestep 



1e-01 



FIG. 5. Error as function of timestep for all timestepping algorithms, for a random medium and random initial conditions, 
using the data from table II. 



random medium (also studied in for example [49]) and starting from a random initial condition. The results are also 
shown in figure 5. From the table and the figure, it is clear that in this case, the LTS family of algorithms perform 
better than the VS-FDTD algorithms. Again we see that the error scales according to the order of accuracy in time, 
except for very small errors (below 10 -10 ). Especially for the LTS-4 algorithm, we see that accumulation of rounding 
errors, due to large number of operations, increases the error. 



VI. CONCLUSIONS 



In this paper, we have introduced algorithms to solve the time-dependent elastodynamic equations, based on a 
matrix exponential approach. The conservation of the underlying skew-symmetry of the first-order partial differential 
equations while discretizing the spatial operators and fields, offers a sound starting point to expand the time evolution 
operator in Chebyshev polynomials. The resulting one-step algorithm is accurate up to machine precision, and this 
statement is justified by rigorous bounds on the error of the unconditionally stable algorithms. The latter class of 
algorithms proved particularly useful if the total energy should be conserved or if a random initial condition is used. 

Finally, the original VS-FDTD algorithm is modified by recasting it into an exponent operator form. In this new 
formulation, the staggered-in-time nature is removed and higher-order in time algorithms are derived, based on the 
lower-order algorithms. Existing VS-FDTD codes can be modified with minor effort to benefit from these advantages. 

In this paper, only free and rigid boundary conditions are considered. Future research is aimed at incorporating 
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absorbing boundary conditions, and the presence of visco-elastic materials. More sophisticated discretization schemes, 
conserving the skew-symmetry, can be easily incorporated [19], and do not require conceptual changes. 
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APPENDIX: DISCRETIZATION AND DECOMPOSITION OF H FOR TWO-DIMENSIONAL 

ISOTROPIC ELASTIC SOLIDS 



In this appendix, it is shown how the matrix 



n 



o 







(41) 



is discretrized conserving its skew-symmetric properties, for two dimensional isotropic elastic solids. In two dimensional 
P-SV wave propagation, no dependency upon y is assumed. For isotropic elastic solids, the stiffness matrix is given 
in terms of the Lame coefficients A and \i by [24] 



A + 2/i A 
C = | A A • 2// 

yU 



(42) 



in the basis a = (a xx ,a zz , 
needs VC, which reads 




Oxz) T ■ In the skew-symmetric basis, using the variables w = y/pv and s = C x ^ 2 o, one 



(43) 



where 

a = 71 (v /A + A* + V7i) , 

This gives for the explicit form of the matrix 7i 
( 



(44) 
(45) 



H = 



_LJL a _LiL« J_iL m 

y/pdx a ^pdx> J JpdzVr 



a dx VP 

3-9- ±- 



(46) 



in the basis ip = (s xx , s zz , s xz , w x ,w z ) T . The discrete analogue of ip, the vector 'J, is obtained by mapping the fields 
onto the two-dimensional staggered grid [8] that is shown in figure 6. We adopt the convention that the s xx and 
s zz stress fields are located in the (1, 1) corner of the grid. The values of the discretized fields / are related to their 
continuous counterparts g by 



f(i,j,t) = g(i6/2,j5/2,t). 
Therefore, using the unit vector e(i,j, k), fields within the vector \I/ can be indexed on the grid by 



(47) 



*(«,J, s xx ,t) = e£. = s xx (i,j,t), 



(48) 
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FIG. 6. Unit cell of the two dimensional staggered grid onto which the continuous velocity and stress fields of the elastody- 
namic equations are mapped in order to conserve the skew-symmetry. 



and analogous equations apply for indexing the other stress and velocity fields within 

Due to the staggered nature of the grid and the choice of the origin, the s xx and s zz stress fields are only defined 
on the x =odd and z =odd lattice points. Similarly, the w x and w z velocity fields are defined at respectively the 
x =even/z =odd and x =odd/z =even lattice points, and the s xz stress field is given at the x =even and z =even 
grid entries. Note that all for simplicity of notation the fields are indexed on the full grid, despite the fact that they 
are not defined on each point. 

It is assumed that the total number of lattice points in each direction is odd, and also that the boundary is located 
at the first and last rows/columns of the grid. The free or rigid boundary conditions themselves are implemented by 
excluding the field points that are located at the boundary and should remain zero during the time integration. 

Using this grid and the central-difference approximation to the spatial derivative, we obtain the spatially discretized 
analogue of equation (46), for example 



9 , ,, .,1 



w x {i + l,j,t) w x (i-l,j,t) 



and 



d_ 

dt 



w x (i,j,t) = 



a(i + l,j)s xx (i + l,j,t) - a(i - l,j)s xx (i - l,j,t) 
(3(i + l,j)s zz (i + l,j,t) - P(i - l,j)s zz (i - l,j,t) 



(49) 



VW + l)s X z(i,j + !,<)- ~ l)(s xz (i,j - l,t) 

Similar equations hold for s zz , s xz and w z . Using this notation, the matrix H can be decomposed into 

Y[ = Jj(x,S xx ,W x ) _|_ Jj(x,Szz,W x ) _|_ ff(x,S xz ,Wz) _|_ Jj{z,S xx ,W x ) _|_ Jj{z,Sz, z ,Wz) _|_ jj(z 

where, for instance, the explicit form of H^ x ' Sxx ' Wx ^ is given by 

-2 



lz.s xz ,w 



>) 



(50) 



(51) 



_^i,j,s xx ^i+l,j,w x ^i-\-X,j,w x ^i,j,s x 



3=1 i= 



a{i,j) 



=1 i=2 Sy/p{l + l,j) 



j,s xx i+l,j,w x 



e i+l,j,w x e i,j,s x 



H 



(x,s xx ,w x ) 



(52) 



(53) 



Here, the prime in the summation indicates that the summation index is increased with strides of two. It is easy 



to convince oneself that the matrices H\ 



and H. 



(x,s xx ,w x ) 



are block diagonal and skew-symmetric. The other 
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matrices in equation (51) have similar explicit forms and can also be decomposed into block diagonal parts. Therefore, 
a first-order approximation to the matrix exponent, as given by equation (17), reads 

exp(rH) = 

exp (r//;"^ ) exp (r//./ ~ ) exp (r//, " " ) exp (tH^'""' w ' ) ) 
•exp {rH[ x '"" v " ) ') exp (rll, ' ) exp (rH[ z ' s ^ w ^ exp (tH^""""^ 
•exp (rH[ z ' s ^ w ^ exp (tJ# , "" , "* ) ) exp (r//,' " " ) exp (r/// " " ) 

+G(r 2 ). (54) 
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